################################################################################
## 
## Main script for reproducing Input Data
## 
## World Development submission:
##
## Estimating the Spillover Economic Effects of Foreign Conflict: 
##      Evidence from Boko Haram
##
## initial 2020
## revised February 2025
##
## Authors: R. Jedwab, B. Blankespoor, T. Masaki and C. Rodr�guez-Castel�n
##
## Code: B. Blankespoor
##
################################################################################

get_current_file_path <- function(){
  this_file <- grep("^--file=", commandArgs(), value = TRUE)
  this_file <- gsub("^--file=", "", this_file)
  if (length(this_file) == 0) this_file <- rstudioapi::getSourceEditorContext()$path
  return(dirname(this_file))
}

# set working directory
wkdir <- get_current_file_path()
print(wkdir)
setwd(wkdir)
getwd()

print(Sys.time())

# libraries, Rstata
source(file.path(wkdir,'bh_spillover_wd_global_libraries.R'))

## Turn off sf's spherical ("S2") backend
sf_use_s2(FALSE) 
today <- Sys.Date()

# Run date
run_vers <- '2025-02-03'

# renv
#renv::init()
tbldir <- file.path(wkdir,'tables_rp',sprintf("vers%s",run_vers)) # reproduced tables
figdir <- file.path(wkdir,'figures_rp',sprintf("vers%s",run_vers)) # reproduced figures
datdir <- file.path(wkdir,'local') # data directory
prcdir <- file.path(wkdir,'proc_data') # processed data directory
xfun::dir_create(tbldir)
xfun::dir_create(figdir)
xfun::dir_create(datdir)
xfun::dir_create(prcdir)

print("Requires Stata, Python with arcpy and R")
print("User defined directories required")

# stata path here #
stata_path = '"C:/Program Files/Stata18/StataMP-64.exe"'

# python path with arcpy library from ESRI ArcGIS Pro
python_path = "C:/Program Files/ArcGIS/Pro/bin/Python/envs/arcgispro-py3/python.exe"

# conda path
conda_path <- "/AppData/Local/r-miniconda/Scripts/conda.exe"

##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "//Rtools//mingw_$(WIN)//bin//")

## GRID identifier : GRID.dta
GRID_fn <- file.path(wkdir,'proc_data','GRID.dta')
if(!file.exists(GRID_fn)) {
  py_file_path = file.path(wkdir,'fishnet_xy_intersection.py')
  system2('conda activate arcgispro-py3') # activate
  system(paste(python_path, py_file_path),show.output.on.console = TRUE)
  do_file_path = file.path(wkdir,'do_prep_bh_grid.do')
  system(paste(stata_path, "/e do", do_file_path),show.output.on.console = TRUE)
}

## Treatment area dummies
TREAT_fn <- file.path(wkdir,'proc_data','TREAT.dta')
if(!file.exists(TREAT_fn)){
  do_prep_grid_treatment_fn <- file.path(wkdir,'bh_spillover_do_prep_grid_treatment.do')
  system(paste(stata_path, "/e do", do_prep_grid_treatment_fn),show.output.on.console = TRUE)
}

## Niger communes : NER_ADM4.dta
NER_ADM4_fn <- file.path(wkdir,'proc_data','NER_ADM4.dta')
if(!file.exists(TREAT_fn)){
  do_prep_ner_adm4_fn <- file.path(wkdir,'do_prep_ner_adm4.do')
  system(paste(stata_path, "/e do", do_prep_ner_adm4_fn),show.output.on.console = TRUE)
}


##
## lights and pop
##

# Night Time Lights
NTL_POP_ts_fn <- file.path(wkdir,'proc_data','NTL_POP_ts.dta')
if(!file.exists(NTL_POP_ts_fn)){
  source(file.path(wkdir,'prep_zonal_extract_boko_haram.R'))
  do_prep_bh_ntl_pop_acled_fn <- file.path(wkdir,'do_prep_bh_ntl_pop_acled.do')
  system(paste(stata_path, "/e do", do_prep_bh_ntl_pop_acled_fn),show.output.on.console = TRUE)
}

# Harmonized Night Time Lights
NTLH9218_fn <- file.path(wkdir,'proc_data','NTLH9218.dta')
if(!file.exists(NTLH9218_fn)){
  ntlh_dir = '010_imageryBaseMaps/NTL_harmonized' # directory with NTLH
  source(file.path(wkdir,'bh_spillover_prep_ntl_harmonized.R'))
  do_prep_ntl_harmonized_fn <- file.path(wkdir,'do_prep_ntl_harmonized.do')
  system(paste(stata_path, "/e do", do_prep_ntl_harmonized_fn),show.output.on.console = TRUE)
}

# VIIRS Night Time Lights
NTL_VIIRS_vars_1220_fn <- file.path(wkdir,'proc_data','NTL_VIIRS_vars_1220.dta')
if(!file.exists(NTL_VIIRS_vars_1220_fn)){
  # Preprocessed files provided else download 
  source(file.path(wkdir,'prep_ntl_viirs_server_ts.R'))
  do_prep_viirs_v2_fn <- file.path(wkdir,'do_prep_viirs_v2_from_tifs.do')
  system(paste(stata_path, "/e do", do_prep_viirs_v2_fn),show.output.on.console = TRUE)
}


##
## conflict
##

# ACLED
ACLED_w_types_0018_ts3_fn <- file.path(wkdir,'proc_data','ACLED_w_types_0018_ts3.dta')
if(!file.exists(ACLED_w_types_0018_ts3_fn)){
  source(file.path(wkdir,'prep_acled_data_all_1999_2015.R'))
  do_prep_bh_acled_conflict_w_type_grid_fn <- file.path(wkdir,'do_prep_bh_acled_conflict_w_type_grid.do')
  system(paste(stata_path, "/e do", do_prep_bh_acled_conflict_w_type_grid_fn),show.output.on.console = TRUE)
}

# SCAD
SCAD_loc_level_fn <- file.path(wkdir,'proc_data','SCAD_loc_level.dta')
if(!file.exists(SCAD_loc_level_fn)){
  source(file.path(wkdir,'prep_scad.r'))
  do_prep_bh_scad_fn <- file.path(wkdir,'do_prep_bh_scad.do')
  system(paste(stata_path, "/e do", do_prep_bh_scad_fn),show.output.on.console = TRUE)
}

# UCDP
UCDP_loc_level_fn <- file.path(wkdir,'proc_data','UCDP_loc_level.dta')
if(!file.exists(UCDP_loc_level_fn)){
  source(file.path(wkdir,'prep_conflict_ucdp.r'))
  do_prep_ucdp_fn <- file.path(wkdir,'do_prep_ucdp.do')
  system(paste(stata_path, "/e do", do_prep_ucdp_fn),show.output.on.console = TRUE)
}

##
## distance variables
##

# distance variables 
DIST_fn <- file.path(wkdir,'proc_data','DIST.dta')
if(!file.exists(DIST_fn)){
  ## downloaded .tifs of stable lights product
  ## https://www.ncei.noaa.gov/products/dmsp-operational-linescan-system
  ntldmsp_dir <- '/010_imageryBaseMaps/DMSP_stable_lights'
  afripolis_dir = '/016_society/Africapolis_2015'
  source(file.path(wkdir,'prep_zonal_extract_boko_haram.R')) # 2020-09-02
  do_prep_distvars_fn <- file.path(wkdir,'do_prep_distvars.do')
  system(paste(stata_path, "/e do", do_prep_distvars_fn),show.output.on.console = TRUE)
}

# more distance variables 
DIST_CMR_fn <- file.path(wkdir,'proc_data','DIST_CMR.dta')
if(!file.exists(DIST_CMR_fn)){
  afripolis_dir = '/016_society/Africapolis_2015'
  petro_dir = "/016_society/PETRO_dataset_080907/PETRO_Onshore_080907"
  source(file.path(wkdir,'prep_zonal_extract_boko_haram_dist.R')) # 2020-09-02
  do_prep_dist_cmr_fn <- file.path(wkdir,'do_prep_dist_cmr.do')
  system(paste(stata_path, "/e do", do_prep_dist_cmr_fn),show.output.on.console = TRUE)
}

## ROAD-BASED DISTANCE TO BH AREA BORDER AND TO MAIN CITIES
tt2bhbr_fn <- file.path(wkdir,'proc_data','tt2bhbr.dta')
TT_CONTROL_CITIES_fn <- file.path(wkdir,'proc_data','TT_CONTROL_CITIES.dta')
if(!file.exists(tt2bhbr_fn) | !file.exists(TT_CONTROL_CITIES_fn)) {
  # BH BORDER and MAIN CITIES #
  source(file.path('bh_spillover_wd_travel_time2port.R'))
  do_prep_tt2bhbr_fn <- file.path(wkdir,'do_prep_tt2bhbr.do')
  system(paste(stata_path, "/e do", do_prep_tt2bhbr_fn),show.output.on.console = TRUE)
  do_prep_bh_tt_city_controls_fn <- file.path(wkdir,'do_prep_bh_tt_city_controls.do')
  system(paste(stata_path, "/e do", do_prep_bh_tt_city_controls_fn),show.output.on.console = TRUE)
}


## LAKE CHAD AREA
WBLKC_fn <- file.path(wkdir,'proc_data','WBLKC.dta')
if(!file.exists(WBLKC_fn)){
  source(file.path(wkdir,'prep_lake_chad_boundary.R'))
  do_prep_world_bank_lake_chad_boundary_fn <- file.path(wkdir,'do_prep_world_bank_lake_chad_boundary.do')
  system(paste(stata_path, "/e do", do_prep_world_bank_lake_chad_boundary_fn),show.output.on.console = TRUE)
}

## DISTANCE TO LAKE CHAD
DIST2LAKECHAD_fn <- file.path(wkdir,'proc_data','DIST2LAKECHAD.dta')
if(!file.exists(DIST2LAKECHAD_fn)){
  source(file.path(wkdir,'dist2lakechad.R'))
  do_prep_world_bank_lake_chad_boundary_fn <- file.path(wkdir,'do_prep_world_bank_lake_chad_boundary.do')
  system(paste(stata_path, "/e do", do_prep_world_bank_lake_chad_boundary_fn),show.output.on.console = TRUE)
}

## BURNED AREA
BURN_mo_ts_fn <- file.path(wkdir,'proc_data','BURN_mo_ts.dta')
if(!file.exists(BURN_mo_ts_fn)){
  # requires output from GEE
  do_prep_bh_burn_month_ts_fn <- file.path(wkdir,'do_prep_bh_burn_month_ts.do')
  system(paste(stata_path, "/e do", do_prep_bh_burn_month_ts_fn),show.output.on.console = TRUE)
}

## GREENNESSS
NDVI_mo_ts_fn <- file.path(wkdir,'proc_data','NDVI_mo_ts.dta')
if(!file.exists(NDVI_mo_ts_fn)){
  # requires output from GEE
  do_prep_bh_burn_ndvi_ts_fn <- file.path(wkdir,'do_prep_bh_burn_ndvi_ts.do')
  system(paste(stata_path, "/e do", do_prep_bh_burn_ndvi_ts_fn),show.output.on.console = TRUE)
}

## TYPE OF LAND USE
LULCfreq_ts_v2_fn <- file.path(wkdir,'proc_data','LULCfreq_ts_v2.dta')
if(!file.exists(LULCfreq_ts_v2_fn)){
  # requires download
  lulc_dir = '/007_environment/landcover/ESACCI-LC' # ESA data directory
  source(file.path(wkdir,'prep_bh_landcover_reclass_crosstab_1992_2019.R'))
  do_bh_prep_lulc_fn <- file.path(wkdir,'do_bh_prep_lulc.do')  # VELOX PACKAGE NO LONGER SUPPORTED
  system(paste(stata_path, "/e do", do_bh_prep_lulc_fn),show.output.on.console = TRUE) # VELOX PACKAGE NO LONGER SUPPORTED
}

LULC08_share_fn <- file.path(wkdir,'proc_data','LULC08_share2024-07-19.dta')
if(!file.exists(LULC08_share_fn)){
  lulc_dir = '/007_environment/landcover/ESACCI-LC' # ESA data directory
  source(file.path(wkdir,'bh_spillover_wd_prep_landcover_reclass_2008_only_vector.R'))
}

## BORDER CTRLS
GRID_BORDER_fn <- file.path(wkdir,'proc_data','GRID_BORDER.dta')
if(!file.exists(GRID_BORDER_fn)){
  do_prep_grid_borders_fn <- file.path(wkdir,'do_prep_grid_borders.do')
  system(paste(stata_path, "/e do", do_prep_grid_borders_fn),show.output.on.console = TRUE)
}

## MINES, OILREF, OIL PROD
OILPROD_fn <- file.path(wkdir,'proc_data','OILPROD.dta')
OILREF_fn <- file.path(wkdir,'proc_data','OILREF.dta')
URANMINE_fn <- file.path(wkdir,'proc_data','URANMINE.dta')
if(!file.exists(OILPROD_fn) | !file.exists(OILREF_fn) | !file.exists(URANMINE_fn)){
  do_prep_petro_onshore_in_fishnet_fn <- file.path(wkdir,'do_prep_petro_onshore_in_fishnet.do')
  system(paste(stata_path, "/e do", do_prep_petro_onshore_in_fishnet_fn),show.output.on.console = TRUE)
}
 

## REFUGEES
IDPREF_fn <- file.path(wkdir,'proc_data','IDPREF.dta')
if(!file.exists(IDPREF_fn)){
  do_prep_idpref_in_fishnet_fn <- file.path(wkdir,'do_prep_idpref_in_fishnet.do')
  system(paste(stata_path, "/e do", do_prep_idpref_in_fishnet_fn),show.output.on.console = TRUE)
}


## ROADS IN CELL AND DISTANCE TO ROADS
ROADS_fn <- file.path(wkdir,'proc_data','ROADS.dta')
ROADSTYPEM_fn <- file.path(wkdir,'proc_data','ROADTYPEM_ts02262021.dta')
if(any(!file.exists(c(ROADS_fn,ROADSTYPEM_fn)))){
  prep_road_measures_py = file.path(wkdir,'prep_road_measures.py')
  system2('conda activate arcgispro-py3') # activate
  system(paste(python_path, prep_road_measures_py),show.output.on.console = TRUE)
  
  source(file.path(wkdir,'road_length.R'))
  
  # intersections
  do_prep_road_intersections_fn <- file.path(wkdir,'do_prep_road_intersections.do')
  system(paste(stata_path, "/e do", do_prep_road_intersections_fn),show.output.on.console = TRUE)
  
  # ROAD TYPE STOCKS
  do_prep_road_length_ts_fn <- file.path(wkdir,'do_prep_road_length_ts.do')
  system(paste(stata_path, "/e do", do_prep_road_length_ts_fn),show.output.on.console = TRUE)
}


# DISTANCE TO ROADS
DIST2ROADS_fn <- file.path(wkdir,'proc_data','DIST2ROADS.dta')
if(!file.exists(DIST2ROADS_fn)){
  source(file.path(wkdir,'dist2roads.R'))
  do_prep_dist2road_type_fn <- file.path(wkdir,'do_prep_dist2road_type.do')
  system(paste(stata_path, "/e do", do_prep_dist2road_type_fn),show.output.on.console = TRUE)
}

## AIRPORTS
AIRPORT_fn <- file.path(wkdir,'proc_data','AIRPORT.dta')
if(!file.exists(AIRPORT_fn)){
  source(file.path(wkdir,'do_prep_airports.R'))
  do_prep_airports_fn <- file.path(wkdir,'do_prep_airports.do')
  system(paste(stata_path, "/e do", do_prep_airports_fn),show.output.on.console = TRUE)
}

## COTTON
COTT_fn <- file.path(wkdir,'proc_data','COTT.dta')
if(!file.exists(AIRPORT_fn)){
  source(file.path(wkdir,'cotton_production.R'))
  do_prep_cotton_fn <- file.path(wkdir,'do_prep_cotton.do')
  system(paste(stata_path, "/e do", do_prep_cotton_fn),show.output.on.console = TRUE)
}

## ELECTRICITY
ELEC_all_fn <- file.path(wkdir,'proc_data','ELEC_all.dta')
if(!file.exists(ELEC_all_fn)){
  source(file.path(wkdir,'dist2electricity.R'))
  do_prep_dist2elec_fn <- file.path(wkdir,'do_prep_dist2elec.do')
  system(paste(stata_path, "/e do", do_prep_dist2elec_fn),show.output.on.console = TRUE)
}

## POWERPLANT
DIST2PP_all_fn <- file.path(wkdir,'proc_data','DIST2PP_all.dta')
if(!file.exists(DIST2PP_all_fn)){
  source(file.path(wkdir,'dist2powerplant.R'))
  do_prep_dist2powerplant_fn <- file.path(wkdir,'do_prep_dist2powerplant.do')
  system(paste(stata_path, "/e do", do_prep_dist2powerplant_fn),show.output.on.console = TRUE)
}

## POLITICAL CITIES
POLICITYDIST_fn <- file.path(wkdir,'proc_data','POLICITYDIST.dta')
if(!file.exists(POLICITYDIST_fn)){
  source(file.path(wkdir,'do_prep_political_city.R'))
  do_prep_political_cities_fn <- file.path(wkdir,'do_prep_political_cities.do')
  system(paste(stata_path, "/e do", do_prep_political_cities_fn),show.output.on.console = TRUE)
}


## BUILT POP
BUILT_POP_fn <- file.path(wkdir,'proc_data','BUILT_POP.dta')
if(!file.exists(BUILT_POP_fn)){
  ghs_dir = '/016_society/GHS_BUILT/GHS_BUILT_LDSMT_GLOBE_R2018A'
  ghspop_dir = '/016_society/GHS_POP/POP250_R2019'
  source(file.path(wkdir,'prep_ghs_pop_and_builtup.R'))
  do_prep_ghs_new_fn <- file.path(wkdir,'do_prep_ghs_new.do')
  system(paste(stata_path, "/e do", do_prep_ghs_new_fn),show.output.on.console = TRUE)
}


## ETHNIC AREA - GREG, Murdock
ETHNICgreg_fn <- file.path(wkdir,'proc_data','ETHNICgreg.dta')
ETHNICmurd_fn <- file.path(wkdir,'proc_data','ETHNICmurd.dta')
if(!(file.exists(ETHNICgreg_fn) | file.exists(ETHNICgreg_fn))) {
  prep_ethnic_greg_fn <- file.path(wkdir, 'prep_ethnic_greg.py')
  system(paste(conda_path, "activate arcgispro-py3 &&", python_path, prep_ethnic_greg_fn), intern = TRUE, wait = TRUE)
  system2(python_path, prep_ethnic_greg_fn)
  do_prep_ethnic_v2_fn <- file.path(wkdir,'do_prep_ethnic_v2.do')
  system(paste(stata_path, "/e do", do_prep_ethnic_v2_fn),show.output.on.console = TRUE)
}


## GROUNDNUT PROD : GNUT
GNUT_fn <- file.path(wkdir,'proc_data','GNUT.dta')
if(!file.exists(GNUT_fn)){
  # ground nut data
  gnut_dir <- paste0('/001_farming/SPAM_2010_v2r0')
  source(file.path(wkdir,'groundnut_production.R'))
  do_prep_groundnut_fn <- file.path(wkdir,'do_prep_groundnut.do')
  system(paste(stata_path, "/e do", do_prep_groundnut_fn),show.output.on.console = TRUE)
}

## AG SUITABILITY
AGSUIT_fn <- file.path(wkdir,'proc_data','AGSUIT2020-10-20.dta')
if(!file.exists(AGSUIT_fn)){
  afclima_dir <- '/001_farming/Climaf' # directory
  source(file.path(wkdir,'prep_ag_suit_afclima.R'))
  do_prep_ag_suit_v2_fn <- file.path(wkdir,'do_prep_ag_suit_v2.do')
  system(paste(stata_path, "/e do", do_prep_ag_suit_v2_fn),show.output.on.console = TRUE)
}


## COMBINED SUITABILITY
LANDSUIT_fn <- file.path(wkdir,'proc_data','LANDSUIT_share2024-07-19.dta')
if(!file.exists(LANDSUIT_fn)){
  source(file.path(wkdir,'sh_area_land_suitability.R'))
}

## BORDER CROSSINGS
DIST2BCROSS_fn <- file.path(wkdir,'proc_data','DIST2BCROSS.dta')
if(!file.exists(DIST2BCROSS_fn)){
  source(file.path(wkdir,'dist2border_post_in_country.r'))
  do_prep_dist2border_crossing_fn <- file.path(wkdir,'do_prep_dist2border_crossing.do')
  system(paste(stata_path, "/e do", do_prep_dist2border_crossing_fn),show.output.on.console = TRUE)
}

## CROSSWALK
CROSSWALK_fn <- file.path(wkdir,'proc_data','CROSSWALKsorted.dta')
if(!file.exists(CROSSWALK_fn)){
  source(file.path(wkdir,'prep_crosswalk_grid01_x_adm1.R'))
  do_prep_crosswalk_fn <- file.path(wkdir,'do_prep_crosswalk.do')
  system(paste(stata_path, "/e do", do_prep_crosswalk_fn),show.output.on.console = TRUE)
}
 
## MARKET ACCESS
MA_NTL_TT_HR08_fn <- file.path(wkdir,'proc_data','MA_NTL_TT_HR08.dta')
ma_ntl_excl20km_bhbr_shocks_fn <- file.path(wkdir,'proc_data','ma_ntl_excl20km_bhbr_shocks_newlabels200820240725.dta')
if(any(!file.exists(ma_ntl_excl20km_bhbr_shocks_fn,ma_ntl_excl20km_bhbr_shocks_fn))){
  ntlharmdir = '/010_imageryBaseMaps/NTL_harmonized'
  ntldir = '/010_imageryBaseMaps/DMSP_stable_lights'
  slope_dir <- "/006_elevation/slope/meanslope"
  roadsafr_dir <- '/018_transportation/afrroads'
  source(file.path(wkdir,'prep_ma_ntl_tt_excl20km_server_ts_2021_0319_wp.r'))
}

## GROUNDWATER DEPTH
GWDEPTH_fn <- file.path(wkdir,'proc_data','GWDEPTH_share2024-07-18.dta')
if(!file.exists(GWDEPTH_fn)){
  source(file.path(wkdir,'sh_area_gw_depth_cat.R'))
}

## LIVESTOCK PRODUCTION ##
BUFFALO_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_buffalo_2024-07-21.dta')
CATTLE_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_cattle_2024-07-21.dta')
CHICKENS_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_chickens_2024-07-21.dta')
GOATS_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_goats_2024-07-21.dta')
HORSES_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_horses_2024-07-21.dta')
PIGS_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_pigs_2024-07-21.dta')
SHEEP_TLU_fn <- file.path(wkdir,'proc_data','LIVESTOCK_TLU_sheep_2024-07-21.dta')

if(any(!file.exists(c(BUFFALO_TLU_fn,CATTLE_TLU_fn,CHICKENS_TLU_fn)))){
   source(file.path(wkdir,'livestock_production.R'))
}

## ** TRANSCONTINENTAL ROADS **
DIST2TRANS_IMP_fn <- file.path(wkdir,'proc_data','DIST2TRANS_IMP_roads_2024-07-22.dta')
if(!file.exists(DIST2TRANS_IMP_fn)){
  source(file.path(wkdir,'dist2transcontinental_roads.R'))
}

## WOODY BIOMASS **
WOODYBIOMASS_fn <- file.path(wkdir,'proc_data','WOODYBIOMASS_2024-07-22.dta')
if(!file.exists(WOODYBIOMASS_fn)){
  woody_biomass_dir <- '/007_environment/woody_biomass'
  source(file.path(wkdir,'woody_biomass.R'))
}

##
## Figure 1: Boko Haram Area and the Three Countries of Study
##

# updated figure with 3 BH states
source(file.path(wkdir,'bh_spillover_cumulative_events_states.R'))

##
## Figure 2:
##
source(file.path(wkdir,'prep_ACLED_source_figures_until2018_bw.R'))
